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The formulation of a one-dimensional analytical model 
representing steady state operation of a vertically oriented 
gas loaded variable conductance heat pipe is presented. The 
model includes not only the effects of binary mass diffusion 
and axial pipe-wall conduction, but also gravitational 
effects which exist at the vapor-gas interface region when 
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integration using a variable step fourth order Runge Kutta 
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INTRODUCTION 



I . 



A. BACKGROUND 

The heat pipe is a heat transfer device capable of trans- 
ferring large quantities of heat under nearly isothermal 
operating conditions. This device is basically a liquid-vapor 
two phase system which transfers heat by evaporation and 
condensation of a working fluid at opposing ends of a tube. 
There are no moving parts in a heat pipe. To sustain opera- 
tion the condensed working fluid is transferred from the 
condenser to the evaporator section of the device by capillary 
pumping action through a wick structure lining the heat pipe 
inner wall. 

B. CONVENTIONAL HEAT PIPE PRINCIPLES 

The conventional heat pipe in its simplest form consists 
of a hollow tube with wicking material lining its inner sur- 
face as shown in Figure 1. This device is evacuated, then 
filled with an appropriate amount of a working fluid. 

In operation, heat at a given rate from a source is ab- 
sorbed at one end of the pipe evaporating the fluid from the 
wick. The vapor then flows along the pipe to the condenser 
end where heat is rejected to a sink and the vapor is condensed 
into the wick. Circulation of the working fluid is maintained 
by a capillary head in the wick which must exactly balance 
pressure drops in the liquid and vapor, as well as any body 
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Figure 1. Conventional Heat Pipe 
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forces present. The heat load alone controls the satura- 
tion temperature and pressure of the working fluid, which 
become the operating temperature and pressure of the pipe. 

Since the transfer of heat is effected by adding and 
removing the heat of vaporization at saturation conditions, 
the heat pipe is nearly isothermal along its length. How- 
ever, a change in the heat load will automatically change 
the operating temperature (and pressure) of the pipe, given 
a constant sink temperature. This characteristic is unde- 
sirable for many applications where a constant source tem- 
perature is required. 

C. VARIABLE CONDUCTANCE HEAT PIPE PRINCIPLES 

1 . General 

If a noncondensible gas is introduced into a heat 
pipe along with the working fluid, the heat pipe becomes a 
passively controlled device which can operate at various heat 
loads while maintaining a more nearly constant source 
temperature . 

During operation the working fluid behaves as described 
above; however, the noncondensible gas is entrained in the 
vapor flow and forms a gas plug at the condenser end of the 
pipe. This plug acts as a diffusion barrier to the flowing 
vapor so that the effective condenser length, and therefore 
condenser heat transfer area, is reduced. For a given heat 
load and noncondensible gas loading, the pipe will operate 
at a specific temperature. If the heat load is increased the 
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saturation temperature and partial pressure of the working 
fluid, as well as total pressure in the pipe, tend to in- 
crease. In response, the gas plug is compressed which 
increases the effective condenser heat transfer area, thus 
reducing the system temperature rise. A thorough presentation 
on the theory and design of variable conductance heat pipes 
has been provided by Marcus [1]. 

2 . Flat Front Theory 

Bienert [2] developed the earliest model used to 
examine the operating characteristics of a gas loaded heat 
pipe. This flat front theory assumed that the noncondensible 
gas was in the form of a plug with a sharply defined vapor- 
gas interface perpendicular to the axis of the pipe. This 
simple model neglected axial conduction in the heat pipe wall 
and the diffusion of vapor into the stagnant gas plug. This 
early model provided much insight into the operating charac- 
teristics of gas loaded heat pipes. However, experimental 
results [1] have* shown the model to be inaccurate, primarily 
because it neglects axial conduction and diffusion. 

3 . Diffuse Front Theory 

A more accurate model of the gas loaded heat pipe was 
provided by Edwards and Marcus [3]. The model, referred to 
as the 'diffuse front theory', includes effects of binary 
diffusion between the working fluid and the noncondensible gas 
as well as axial conduction in the pipe wall. Marcus reported 
[1] good correlation between experimental results and a com- 
puter program based upon the diffuse front theory. 
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4. Gravitational Effects 



Neither the flat front theory nor the diffuse .front 
theory included body force effects on the vapor-gas inter- 
face region. This simplification was probably made because 
much of the early work with heat pipes was related to space 
program applications. 

Kelleher and Batts [4] reported convincing experimen- 
tal results indicating that gravity effects can significantly 
distort the axial temperature profile when the working fluid 
and noncondensible gas are of significantly different molecular 
weights . 

D. OBJECTIVE 

The objective of this research was to develop a one-dimen- 
sional analytical model describing steady state operation of 
a gas loaded variable conductance heat pipe. The model was 
to include not only the effects of axial conduction and binary 
mass diffusion, but also gravitational effects on the vapor- 
gas interface region. An additional purpose of this effort 
was to transform the analytical equations into a set of non- 
dimensionalized simultaneous differential equations which 
could readily be solved numerically using IBM's Continuous 
System Modeling Program III [5]. 
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II . FORMULATION OF ANALYTIC EQUATIONS 



In formulating the analytic equations describing steady 
state heat pipe operation the differential control volume 
in Figure 2 was developed to represent a section of the heat 
pipe condenser. As implied in Figure 2, only vertically 



Q 




q + dq 



p a v a + dS? ( Pa v a )dx 



Figure 2. Differential Volume Element 
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oriented heat pipes are considered so that the transport 
processes can be treated as one dimensional. 

A. SIMPLIFYING ASSUMPTIONS 

The following assumptions were made in formulating an 
analytical model using the various transport processes: 

1. The heat, mass, and momentum transport processes were 
assumed to be one-dimensional in the axial direction within 
the vapor space. This assumption restricted this analytical 
model to the specific case of a vertically oriented heat pipe. 
In any other orientation, gravitational effects should be 
expected to cause highly three dimensional spatial variations 
when the molecular weights of the working fluid and noncon- 
densible gas differ substantially. This effect was demon- 
strated experimentally by Kelleher and Batts [4]. 

2. The working fluid and the noncondensible gas were treated 
as thermodynamically independent substances. 

3. Steady state conditions were assumed. 

4. The working fluid was considered to be at saturation 
conditions at all times. 

5. Thermal resistance of the wick-wall structure, in the 
radial direction, was considered to be negligible compared 
to the thermal resistance outside of the pipe. 

6. The noncondensible gas was treated as a perfect gas. 
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B. CONSERVATION OF MASS 



1 . Species 'a' (working fluid) 

Referring to the differential volume element in 
Figure 2 and considering the system on a unit time basis, 

mass of 'a' in = mass of 'a' out 



or. 



P 



v Ac 
a a 



Ip a v a 




V 

a a 



) dx] Ac + m Pvdx 
a 



which reduces to 



(II. 1) 



0 



m Pv + Ac " 3 — ( p v ) 
a dx *a a 



(II. 2) 



Fick's first law of binary diffusion, as developed 

in reference 6, was used to relate the mass flux of working 

fluid (p v ) to the mass flux of the mixture (pv) . For 
a a 

the one-dimensional case the scalar form of Fick's law can 
be written. 



p a v a * “a (pv) 



-pD 



ab 




(II. 3) 



where the mass fraction of species 'a' is defined as. 



“a = p a /p 



(II. 4) 



and the density of the mixture is given by 
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p 



(II. 5) 



= p + p, 

w a M b 



since species 'a' and species ' b * occupy the same volume as 
independent substances. 

Then equation II. 2 becomes, 



. Pv , d , . d , dtJ a 

m - — + -= — ( co pv) - (pD , - 3 — 

a Ac dx a dx ab dx 



) = 0 



(II. 6 ) 



which represents the conservation of mass for species 'a'. 

2 . Species 'b 1 (noncondensible gas) 

A similar equation was written for species 'b' noting 
that since species 'b' is a noncondensible gas, m^ E 0 . 

Then , 



dx (w b pv) " dx (pD ba dx ) 0 



(II. 7) 



where the mass fraction of species 'b' is defined as. 



“b = p b /p 



(II. 8 ) 



3 . Mass Continuity of the Mixture 

Adding the equations for species 'a' and species 'b' 
(III . 6 and III. 7) gave the continuity equation for the 
mixture , 



p . . dw dw, 

*.II + ax I< ’ v ‘“a + V 1 -E IpD( ^ + ^ ,! = 0 



dx dx 



(II. 9) 



where D = is used for simplicity of notation. 
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Recalling equations II. 4, II. 5, and II. 8, 



which reduced equation II. 9 for mass continuity of the mix- 
ture to. 



m 



Eiz: 

Ac 



+ fj(PV) 



= 0 



(II. 10) 



C. CONSERVATION OF MOMENTUM 

1 . Mixture of Species 'a' and Species 'b 1 

Conservation of linear momentum for the mixture, with 
steady state conditions, was written for the vapor space 
differential element of Figure 3 as, 

£f = / v(pv • dA) 

cs 

— y 

where, dA is the outward directed surface area element and, 

— y 

F is the sum of all external forces acting on the fluid in 

-y — y 

the control volume. For the one-dimensional case F = F . 

x 

The sum of the body forces and surface forces was seen from 
Figure 3 to be, 



£f = p g Ac dx - t Fv dx - dp Ac (II. 11) 

X 

and the surface integral, representing the net efflux of 
linear momentum in the axial direction from the control 
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PvAc 



dA 

A 



pAc 



A 



V 




x 



N 



TPvdx 



pgAcdx 



m Pvdx 
a 






V 






'if 

[pv + d ( pv) ] Ac 




x + dx 




( p + dp) Ac 



Figure 3. Vapor Space Differential Element 
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volume, was expanded to give: 



2 

J = - pv Ac + ( v + dv) ( pv + d( pv) ) Ac + 



CS , X 



v m Pvdx 
a a . 



= Acpvdv + Acvd ( pv) + v m Pvdx 

3 . 3 



( 11 . 12 ) 



equating equations II. 11 and 11.12, and rearranging yielded. 



pg 



Pv _ dp 
Ac dx 



pv 



dv 

dx 



+ 



d( pv) 
dx 



. Pv 
+ v m - — 
a a Ac 



(11.13) 



For steady state conditions the noncondensible gas remains 
stationary with respect to the coordinates fixed to the heat 
pipe. Therefore the net flux of species 'b ' is identically 
zero (i.e., p^v^ = 0 since v^ = 0) . Noting that 

pV = p a V a + P b V b ' and v b " 0 ' 



pv = p V 
3 3 



or 



v 

3 




Repl3cing v in equstion 11.13 and rearrsnging resulted in 

3 

the momentum equation for the mixture. 



dp 

dx 



d , 2. Pv m a 

pg - -=— ( pv ) - - — v — 
^ dx Ac u> 

3 




(11.14) 
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where the viscous shear force x is some function of velocity. 



2 . Species 'b* 

Since the velocity of species ’b' is zero, the con- 
servation of momentum equation for species 'b' could be 
written directly by inspection as 



in form to equation 11.14 could be written. However the addi- 
tional equation would not be useful since it would not be 
linearly independent of the two momentum equations already 
presented. 

D. ENERGY BALANCE FOR PIPE WALL 

In developing a heat pipe energy balance the wick and 
pipe wall were considered together as a composite. Recall 
that thermal resistance of the wick-wall structure, in the 
radial direction, was assumed to be negligible. Figure 4 
then represents an appropriate differential element upon 
which to develop the energy balance. 

By the conservation of energy (referring to Figure 4) , 



8 u 

For example, x = — v, if Poiseuille flow is assumed 




(11.15) 



dx 




A separate momentum equation for species 'a', similar 



in 



Q 



out 



then 
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Qc 




Figure 4. Differential Element for Energy Balance 
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Qc + A h- Pv dx = Qc + dQc + h Ps (T - T ) dx 

a t CT co 



a 



(11.16) 



but, by Fourier's law of conduction 



Qc 




(11.17) 



Incorporating equation 11.17 into equation 11.16 and 
rearranging yields the energy balance equation. 



(a) the coefficient of thermal conductivity (k) 
must be selected to represent the composite 
wick-pipe wall structure, and 

(b) the heat transfer coefficient, h, is some 
function of temperature. The correlation 
selected for h(T) depends upon the external 
cooling conditions. 



E. EQUATIONS OF STATE 
1 . Species 'a' 

Since saturation conditions were assumed to exist 
in the working fluid at all times, tabulated saturation 
properties of the working fluid can be used to generate 
functions for partial pressure, partial density and heat of 




0 



(11.18) 



where 
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vaporization of the working fluid as functions of local 
temperature. Then, 



p a = Fa (T) 

cl 



( 11 . 19 ) 



p a = Sa(T) 

CL 



( 11 . 20 ) 



The functions Fa (T) , Sa(T) , as well as h^^(T) repre- 
sent actual tabulated saturation properties of the working 
fluid which can be readily used in a computerized interpo- 
lation scheme. Alternatively Fa (T) , Sa(T), and h^g(T) could 
be approximated by polynomials or other anlaytic functions 
selected to fit the working fluid saturation curves. 

2 . Species 'b' 

As a perfect gas. 



p b = 




3 . Mixture Properties 

For independent substances. 



( 11 . 21 ) 



p = p a + p b 



( 11 . 22 ) 



and. 



P 



P a + 



( 11 . 23 ) 
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Although equation 11.22 is valid based upon Dalton's law 
of partial pressures, it is included here for completeness 
only since it is not linearly independent of equations 
11.14 and 11.15 developed from momentum considerations. 

F. SUMMARY OF EQUATIONS 

The following set of ten independent equations in ten 
unknowns describes steady state operation of a vertically 
oriented gas loaded heat pipe. 




0 



(II. 6) 




0 



(II .10) 





(11.14) 




(11.15) 




0 



(11.18) 



R 



T 



( 11 . 21 ) 



P b M, 



b 



03 



a 



P a /P 

a 



(II. 4) 



P 




(11.23) 
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Fa (T) 



( 11 . 19 ) 



P 



a 



P a = Sa (T) 

a. 



( 11 . 20 ) 



The ten dependent variables (m , v, p, p , p, , p, p , p, , 

cl 3 D 3. D 

w , T) are all functions of axial position. 

3 
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III. TRANSFORMATION TO COMPUTER MODEL 



A. ELIMINATION OF VARIABLES 

The ultimate goal was to obtain two nondimensionalized 
differential equations with nondimensional temperature and 
mass fraction of species 'a' as the sole dependent variables. 

By reducing the number of equations and dependent variables 
analytically the number of computations required for each 
numerical integration step is significantly reduced. This 
approach was taken to enhance the numerical accuracy and 
stability of computer results for the model. With this 
objective in mind, the set of ten equations in ten unknowns 
summarized in Section II. F was reduced to two equations in 
two unknowns by the following steps. 

1. Total pressure (p) appeared in equation 11.14 only and 
therefore could not be eliminated as a dependent variable by 
back substitution. Since total pressure as a function of 
axial position was not of particular interest, p was eliminated 
by assuming that the axial pressure gradient of total pressure 
was negligible (dp/dx = 0) . Since the axial pressure gradient 
of total pressure is primarily created by the hydrostatic 
head of the vapor-gas mixture, this common assumption is 
reasonable. This assumption means that total pressure is 
assumed to remain essentially constant throughout the vapor- 
gas space of the pipe for a given operating power, gas load, 
operating temperature and sink temperature. 
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It is significant to note that although the axial gradi- 



ent of total pressure is negligible, axial variations in 
the partial pressures of both working fluid vapor and noncon- 
densible gas are significant. For example, using methanol 
as the working fluid at typical operating conditions results 
in an axial partial pressure variation of several atmospheres. 

This assumption eliminated equation 11.14 and the variable 
P- 

2. The following simple substitutions were then sequen- 
tially made in the remaining equations. 

a. Equation 11.20 was used to eliminate p . 

cl 



P 



a 



Sa (T) 



b. Equation 11.19 was used to eliminate p . 

cl 



P 



a 



Fa (T) 



c. Equation 11.21 was used to eliminate p^. 



R 



T 



p b M, 



b 



d. Equation 11.23 was used to eliminate p^. 



P 




e. Equation II. 4 was used to eliminate p. 



P 



a _ Sa (T) 



P 



to 



to 



a 



a 
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3. The surviving four equations in four unknowns ( 
v,T,oj ) were: 

cl 



* + §rr(Sa v) 



a Ac 



dx 



d , Sa D 
dx v oj 

a 



doj 




0 



m 

a 



Pv d , Sa v 

Ac dx oj 

a 



0 




_1_ 

OJ 

a 



1)”] 

V 




i) g 




Ps h 
k Ap 



(T - T ) 

oo 



+ 



h Pv 

-12 

k Ap 



m 

a 



0 



4. Equation III. 2 was then used to eliminate in 
equations III.l and III. 4. 



m 



Ac d , Sa v . 
Pv dx oj 



Equation III.l then became: 



d , Sa v> d . d ,Sa D ^ w a, 

dx to dx dx oj dx 

a a 



= 0 



which was integrated directly to yield. 



o /i 1 \ Sa D ^ W a 

Sa V<1 - — » - — -dF 
a a 



= C 



where C is a constant of integration. However, at x = 



(111.1) 

(111. 2) 

(111. 3) 

(111. 4) 
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V 



0 



and 



Sa/co = p 
a o 



then 



C 



dco 



D p 
o o 



dx 



(III. 5) 



where x = 0 was assumed to be at the condenser-reservoir end 

dco ^ 

= 0 is a boundary condition as 
o 



of the pipe. However, 



dx 



will be discussed later. 



Therefore C = 0 



The remaining system of three equations was then. 



Sa v ( 1 - — ) - ^ -r-2. = 0 

to to dx 

a a 



( III. 6) 



5-[Sa( — - l)-£- T] 
dx co^ 



sa(- - l) g 
a 



(III. 3) 



d T Ps h 



dx 



2 k Ap 



(T-T ) - 



h,. Ac , „ 

fg d , Sa v 



k Ap dx co 



I— 1 ) = 0 



(III. 7) 



5. Equation III. 6 was finally solved for velocity and used 
to eliminate v in equation III. 7 to give a system of two 
equations in T and co^ . 



d_ 

dx 






1)9 




Ps h 
k Ap 



(T - T ) 

OO 



h _ Ac 

Sr 

k Ap 



d r D Sa 
dx l to (co -1) 



dco 



dx 



0 



(III. 3) 



(III . 8) 
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To solve these two equations simultaneously it was advan- 
tageous to expand and differentiate equation III. 3 and to 
expand and rearrange equation III. 8 thereby generating a 
system of two second order differential equations in w 
and T. 

After differentiating equation III. 3 a considerable amount 
of algebra led to the following form: 



d T { (— - 1) (T ^ + Sa) } - 



d oj 



, 2 (o 
dx a 



dT 



a^( T_Sa.) 



dx' 



0) 



= ( 



1 - — ) (2 ^ + T 2 + (-^ ~ (T ^ + Sa) 

in dT . 2 dx 2 dx dx 



dT 



0) 



dT 



- 2 2 V { _1_ . j, 

3 dx R U o) l ’ 1 dT dx 

a) a 



dSa dT^ _ Sa 



dto 



0) 



2 dx 



(III. 9) 



Equation III. 8 was manipulated into a similar convenient 
format as: 





Ac h 
^Ap k 



12_ 

(w 



D Sa 

- 1 ) 0 ) 

a 



-} 

a 



Ps h 
k Ap 



(T - T ) 

oo 



Ac h 



fSL 



do) 



Ap k (io -l)w dx 

cl 3. 



dto 

{ - D Sa<-g£)( 



, — + -i) 

00 -1 00 

a a 



+ (Sa dD dSa. dTi 

+ (Sa dT + D dT ) dx > 



(III . 10) 



Equations III. 9 and III. 10 represent the final form of the 
desired equations in a dimensional format. 



32 



B. BOUNDARY CONDITIONS 

To numerically integrate these second order differential 
equations four boundary conditions must be established at 
x = 0. To avoid the obvious singularity when = 1.0 (at 
the evaporator end of the condenser) integration must proceed 
from the condenser-reservoir end toward the evaporator. 
Therefore x = 0 is chosen to be the condenser- reservoir end 
of the heat pipe. 

Appropriate boundary conditions are then as follows. 

1. T = T at x = 0 (III. 11) 

oo 

2. ~ = 0 at x = 0 (III. 12) 

dx 

doi 

3. - 3 -^- =0 at x = 0 (III. 13) 

dx 



4. a} = oj >0 atx=0 (III. 14) 

cl 3 . 

o 

The selection of an initial value of oj is somewhat arbi- 

a 

trary. However, the initial value of to selected, in conjunc- 

cl 

tion with the values of T and T used in the model, indirectly 

0 oo 

fix the amount of noncondensible gas and the operating power 
level of the heat pipe being simulated for a particular 
working fluid - noncondensible gas pair. 

C. NONDIMENS IONALIZ ATI ON OF EQUATIONS 

To nondimensionalize equations III. 9 and III. 10 the 
following dimensionless variables were defined. 
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XS = x/L nondimensional length 

(independent variable) 

T - T 

OO 

Yl = nondimensional temperature 

e co 



where: 



T = T(x) 

T g = prescribed evaporator temperature 
T m = ambient temperature 



Y2 = In to logarithm of the mass fraction of 

a species 'a' 



For ease of notation and computer programming the first 
and second derivatives of the two dependent variables were 
represented as follows. 



Y3 



d (Yl) 
dXS 



DYl 



Y4 



d(Y2) 

dXS 



DY2 



DY3 



d(DYl) 

dXS 



d 2 (Yl) 
dXS 2 



DY4 



d ( DY 2 ) 
dXS 



d 2 ( Y2) 
dXS 2 



Each function of temperature was nondimensionalized with 
respect to the function value at evaporator temperature (T ) . 
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The following nondimensionalized functions of temperature 
were then developed for direct substitution into equations 
III. 9 and III. 10 . 



G1 



D 

D(T ) 
e 



nondimensional binary mass 
diffusion coefficient 



G2 



d(D) 

dT 



< T e - T J 

D(T e ) 



nondimensional temperature 
gradient of the binary mass 
diffusion coefficient 



G3 







nondimensional heat of 
vaporization 



G4 



Sa 

Sa (T ) 
e 



nondimensional working 
fluid density 



G5 = 



d (Sa) (T e" T “ ) 



dT 



Sa (T ) 
e 



G6 = 



d 2 (Sa) (T e ” T <» ) 



dT 



Sa(T e ) 



nondimensional temperature 
gradient of working fluid 
density 



nondimensional curvature of 
working fluid saturation 
density curve 



G7 



h 

h(T ) 
e 



nondimensional convection 
heat transfer coefficient 



Using the nondimensional variables and functions of 
temperature defined above, equations III. 9 and III. 10 were 
transformed into the following nondimensionalized system of 
two equations in two dimensionless dependent variables (Y1 
and Y2) , and their first and second derivatives with respect 
to axial (dimensionless) length. 
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\ 



DY3{ (1 - Exp (Y2) ) ( (KlYl+1) G5 + K1G4 ) } - DY4 ( (KlYl+1) G4 ) 



= ( 1-Exp (Y2) ) {Y3 (Y3 ( 2K1G5+ ( KlYl+1) G6 ) +K2G5) } 



+ Y4 { 2Y3 ( (KlYl+1) G5+K1G4) - G4 ( (KlYl+1) Y4-K2) } (III. 15) 



DY3 (Exp (Y2) -1) - DY4 (K3G1G3G4) = ( Exp ( Y2 ) -1) K4G7Y1 



+ 



G1G3Y4 (K3 (G4 ( 



G2Y3 Y4 Exp ( Y2 ) 
G1 (Exp ( Y2 ) -1 ) 



) + G5Y3) ) 



(III. 16) 



where Kl through K4 are dimensionless constants: 



Kl 



K2 



K3 



K4 




% 9 L 

(— ) 

V R T g 1 

OO 3 £ 



Ac h (T e ) D(T ) Sa (T ) 
. fg e e e 

* k Ap ( T - T ) 
r e 00 



Ps h (T ) 
e 

k Ap 




(III .17) 



In terms of the nondimensional variables, the boundary 
conditions (equations III. 11 through III. 14) became: 



1. Y1 = 0 at XS = 0 (III. 18) 

2. Y2 = Y2 < 0 at XS = 0 (III. 19) 

o 
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3. Y3 = DY1 = 0 at XS = 0 (III. 20) 

4. Y4 = DY2 = 0 at XS = 0 (III. 21) 

where Y2 = In w . As noted previously with the dimensional 
form of this boundary condition (equation III. 14) , the initial 
value of to a selected is arbitrary. However the initial value 
of , in conjunction with the specified values of T g and T^, 
indirectly specifies the quantity of noncondensible gas and 
operating power level of the heat pipe. Since 0 < to <1.0, 
the corresponding value of Y2 must be negative. 
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IV. DISCUSSION AND RECOMMENDATIONS 



A. DISCUSSION 

The Continuous System Modeling Program (CSMP) , developed 
by the IBM Company [5] , is a straightforward computer package 
available to numerically solve equations III. 15 and III. 16. 

The CSMP package includes a number of optional numerical 
integration methods including a variable step fourth-order 
Runge-Kutta method. 

Although equations III. 15 and III. 16 are highly nonlinear 
in the dependent variables Yl and Y2, this system of equations 
can be considered linear in the second derivatives of Yl and 
Y2. Therefore these equations can be written in the following 
form: 



(C0EF1) (DY3) + ( C0EF2 ) (DY4) = FI (IV. 1) 

(C0EF3) (DY3) + ( C0EF4 ) (DY4) = F2 (IV. 2) 



where, 

COEFl = { (l-Exp(Y2) ) ( (K1Y1+1)G5+K1G4) } (IV. 3) 

C0EF2 = -( (K1Y1+1)G4) (IV. 4) 
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FI = ( 1-Exp (Y2) ) { Y 3 ( Y 3 ( 2K1G5+ (K1Y1 + 1 ) G6 ) +K2G5 ) } 

+ Y4 { 2Y3 ( (K1Y1+1) G5+K1G4)-G4 ( (KlYl+l)Y4 - K2) } 

(IV. 5) 

C0EF3 = (Exp ( Y2 ) -1 ) (IV. 6) 

C0EF4 = - (K3G1G3G4) (IV. 7) 



F2 = (Exp (Y2) -1) K4 G7 Yl + G1 G3 Y4 (K3 (G4 



,G2 Y 3 

~ gT 



Y4Exp ( Y2) 
(Exp ( Y2) -1) 



G5Y3) ) 



( IV . 8 ) 



Using these definitions, the nondimensionalized equations 
representing steady state heat pipe operation may be numerically 
integrated using the variable step fourth-order Runge-Kutta 
option of the CSMP routine as follows. 



1 

Yl = / (DYl)dXS 

0 



1 

/ 



0 



d(Yl) 

dXS 



dXS 



(IV. 9) 



1 

Y2 = / ( DY 2 ) dXS 

0 



1 

/ 



0 



d(Y2) 

dXS 



dXS 



(IV. 10) 



1 

Y3 = / ( DY 3 ) dXS 

0 



j 1 si qn dxs 

0 dXS 



(IV. 11) 
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Y4 



dXS 



(IV. 12) 



1 

= / ( DY 4 ) dXS 

0 



1 2 

j d z (Y2 ) 
0 dXS 2 



where , 



DYl = Y 3 



(IV. 13) 



DY2 = Y4 



(IV. 14) 



and by Cramer's rule. 



DY3 



(FI C0EF4 - F2 C0EF2 ) 

( COEFl COEF4 - COEF2 COEF3 ) 



(IV. 15) 



DY 4 



(F2 COEFl - FI COEF 3) 
(COEFl COEF 4 - COEF2 COEF3 ) 



(IV. 16) 



This model can be used to demonstrate gravitational effects 
on the vapor-gas interface region by providing axial profiles 
of temperature and the mass concentration of the working 
fluid. When the molecular weight of the noncondensible gas 
is less than that of the working fluid, the axial profiles 
will shew a fairly small interface region between vapor and 
gas. However, when the molecular weight of the noncondensible 
gas is substantially greater than that of the working fluid, 
the profiles will show a distorted (spread out) interface 
region. 

Qualitative results can be obtained from this model by 
using additional simplifying assumptions and built-in CSMP 
functions as follows: 
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1. Nominal constant values can be selected to represent 
the convection heat transfer coefficient and the binary mass 



diffusion coefficient. For example, 



h (T) 



2 

1.0 BTU/HR/FT for free convection 



D(T) 



2 

0.4 FT /HR for Krypton-Methanol 



or 



D(T) 



2 

1.9 FT /HR for Helium-Methanol 



2. To generate values of working fluid density (Sa) and 



generator can be used to provide quadratic interpolation of 
tabulated saturation properties. Tabulated thermodynamic 
properties of methanol are available in reference 7. The 
CSMP function "SLOPE" can be used to generate values of 



The temperature profile generated by this model for a 
given vertical heat pipe using a particular working fluid- 
noncondensible gas pair, with given input values of T g , 

T , and oj (0) , represents heat pipe operation at a certain 

00 3 . 

power level and total pressure with a specific load of non- 
condensible gas. The power level, total pressure, and mass 
of noncondensible gas represented could be computed as 
follows . 




2 2 

d(Sa)/dT and d (Sa)/dT as functions of temperature. 
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(IV. 17) 



Mass of Gas (lb ) 



m 




1 

Power Level (BTU/HR) = Pv L f m h £ dXS 

0 a f 9 



(IV. 18) 



Total Pressure (psia) 



Sa(TJ 

Fa ( T ) + TT^-d-w (0)) 

<» w (0) a 



R T 



M, 



b 



00 



(IV. 19) 



B . RECOMMENDATIONS 

Based upon preliminary work using CSMP, numerical results 
from this model could be improved by the following measures. 

1. Analytical expressions are recommended in lieu of the 
special CSMP functions to represent working fluid saturation 
properties as functions of temperature. Expressions for the 
first and second derivatives of density as functions of tem- 
perature could then be found analytically and used directly. 
The resulting expressions should yield more accurate values, 
particularly for the second derivative of density, than are 
achievable through multiple use of the CSMP "SLOPE" function 
with tabulated saturation values. 

2. Since the binary mass diffusion coefficient is a fairly 
strong function of temperature, quantitative results in the 
interface region could be improved considerably by use of 

an empirical expression for D(T) instead of a nominal constant 
value. Reference 8 contains appropriate empirical expressions 
for binary mass diffusion coefficients. 
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3. Quantitative results could also be improved somewhat 
by using an appropriate empirical expression for the convec- 
tion heat transfer coefficient as a function of temperature 
instead of a nominal constant value. Reference 9 provides 

a variety of expressions for h(T) depending upon the convec- 
tive condition being simulated. 

4. With the improvements noted above, this model could be 
used in an iterative loop to obtain preliminary quantitative 
design information concerning vertically oriented terrestrial 
heat pipes when the molecular weights of the operating sub- 
stances are substantially different. 

For example, if the working fluid- noncondensible gas pair 
is specified, along with the required heat source (evaporator) 
temperature range, ambient conditions, and the maximum safe 
total pressure for structural integrity of the pipe, then 
equation IV. 19 could be used to determine the value of w (0) 
to use for boundary condition III. 19. Equation IV. 17 could 
then be used to numerically compute the required mass of non- 
condensible gas to load depending upon the heat pipe dimen- 
sions. Equation IV. 18 would then yield the maximum heat 
rejection rate obtainable for these conditions. If this heat 
rejection rate is not adequate for the design requirements, 
the designer could vary any of a number of parameters, such 
as operating temperature, pipe dimensions, total pressure 
allowed, substances used, etc., until an adequate design is 
found by iteration. 
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